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1  Introduction 


Buoyancy  generated  turbulence  plays  an  important  role  in  the  dynamics  of  atmospheric  boundary  layers 
and  engineering  applications  such  as  space  heating  and  cooling,  electronic  equipment  cooling  and  so¬ 
lar  collectors.  While  shear-generated  turbulence  (mechanical  turbulence)  has  been  studied  extensively 
leading  to  the  development  of  reasonably  reliable  engineering  tools  (turbulence  models),  buoyancy¬ 
generated  turbulence  is  relatively  less  scrutinized  [1] .  Many  of  the  engineering  calculations  of  thermally 
dominated  flows  typically  employ  empirical  correlations  involving  non-dimensional  parameters  such  as 
Nusselt,  Grashoff  and  Rayleigh  numbers.  These  empirical  correlations  lack  generality  and  are  usually 
useful  only  in  situations  for  which  they  have  been  calibrated.  Therefore,  in  buoyant  turbulent  flows, 
there  is  certainly  a  need  for  turbulence  models  that  are  based  on  governing  equations  of  the  flow  and 
hence  possess  a  reasonable  degree  of  generality.  Accurate  modeling  of  buoyancy  generated  turbulence 
requires,  at  least,  a  working  knowledge  and  understanding  of  that  process.  Data  gathered  from  labora¬ 
tory  experiments  are  usually  limited  to  integral  quantities  such  as  heat  and  mass  transfer  coefficients. 
While  these  quantities  are  very  useful  for  ultimately  validating  the  models,  they  lack  sufficient  detail 
to  shed  light  on  the  physics  of  buoyant  turbulence.  Direct  numerical  simulation  (DNS),  a  calculation 
wherein  all  of  the  time  and  length  scales  of  turbulence  are  resolved,  can  provide  data  with  the  required 
degree  of  detail  to  study  the  physics  of  buoyancy  generated  turbulence.  Although  DNS  is  usually  re¬ 
stricted  to  simple  geometry  and  modest  Reynolds  numbers,  it  captures  important  universal  aspects  of 
the  physics  and  has  led  to  the  development  of  some  reasonable  models  of  shear-generated  turbulence. 
It  is  only  recently  that  DNS  data  has  been  used  in  the  development  of  buoyant  turbulence  models.  Dol, 
Hanjalic  and  Kenjeres  [2]  and  Girimaji  and  Hanjalic  [3]  have  used  DNS  turbulent  natural  convection 
data  of  Boudjemadi  et  al.  [4]  and  Versteegh  et  al.  [5]  to  study  certain  modeling  assumptions  and 
ultimately  validate  the  models.  The  turbulent  convection  in  this  case  is  between  two  parallel  vertical 
plates  at  different  temperatures  and  this  configuration  gives  rise  to  shear  as  well  as  buoyancy  generated 
turbulence. 
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Our  objective  in  this  paper  is,  using  DNS  data,  to  study  the  Reynolds  stress  and  turbulent  thermal 
flux  in  a  flow  where  buoyancy  is  the  only  source  of  turbulence  production.  The  turbulent  flow  in 
the  present  case  is  the  Rayleigh-Bernard  convection  between  two  horizontal  parallel  plates  at  different 
temperatures.  The  lower  plate  is  hotter  than  the  upper  plate  leading  to  unstable  stratification.  The 
turbulence  is  statistically  homogeneous  along  horizontal  planes.  The  mean  velocity  is  zero  everywhere 
and,  hence,  there  is  no  shear  production  of  turbulence.  Rayleigh-Bernard  convection  offers  a  unique  flow 
situation  in  which  buoyancy-generated  turbulence  can  be  studied  in  isolation,  free  from  the  complicating 
influence  of  shear-generated  turbulence.  The  four  specific  objectives  of  this  paper  are  to:  (i)  to  evaluate 
the  pressure-strain  and  pressure-temperature  gradient  correlation  models  using  DNS  data;  (ii)  to  study 
the  turbulent  transport  of  Reynolds  stress  and  thermal  flux;  (iii)  to  examine  the  various  modeling 
assumptions  made  to  derive  algebraic  stress  models;  and,  (iv)  to  develop  fully  explicit  algebraic  models 
for  Reynolds  stress  and  thermal  flux. 

The  paper  is  organized  as  follows.  In  Section  2,  the  governing  equations  of  Reynolds  stress  and 
turbulent  thermal  flux  are  presented  along  with  their  standard  closure  models.  Algebraic  models  for 
the  Reynolds  stress  and  thermal  flux  are  also  developed.  Description  of  the  numerical  simulation  of 
turbulent  Rayleigh-Bernard  convection  is  given  in  Section  3.  We  use  the  simulation  data  to  evaluate  and 
develop  various  closure  models  in  the  Reynolds  stress  and  turbulent  thermal  flux  evolution  equations  in 
Section  4.  Also  in  Section  4,  we  examine  the  validity  of  the  various  assumptions  made  in  the  algebraic 
stress  methodology.  We  conclude  in  Section  5  with  a  brief  discussion. 

2  Governing  equations  and  model  development 

The  buoyant  turbulent  flow  considered  here  is  governed  by  the  usual  conservation  of  mass  (continuity), 
momentum  (Navier-Stokes)  and  temperature  (energy)  equations  subject  to  the  Boussinesq  approxima¬ 
tion.  The  flow  variables  are  decomposed  into  their  mean  (upper  case  symbols)  and  fluctuating  (lower 
case  symbols)  parts  and  the  equations  are  Reynolds  averaged  (averaging  is  indicated  by  angular  brack¬ 
ets).  The  resulting  Reynolds-averaged  Navier  Stokes  equations  are  given  in  most  standard  text  books 
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and  are  not  presented  here.  These  averaged  equations  contain  two  unclosed  terms,  the  Reynolds  stress 
(( UiUj ))  and  the  turbulent  thermal  flux  (( urf )). 

2.1  Transport  equations  for  (UiUj)  and  ( urf }. 

The  transport  equation  for  Reynolds  stress  in  a  buoyant  turbulent  flow  (with  no  mean  velocity  gradi¬ 
ents)  is  given  by  (Peeters  and  Henkes  [6]) 

+  Uk{uiUj)tk  —  Gij  -  Eij  +  4>ij  +  T>ij,  (1) 

where  the  subscript ,  k  indicates  differentiation  in  the  k  direction.  The  terms,  respectively,  are  the  time 
rate  of  change,  advection,  buoyant  production  (G^),  dissipation  pressure-strain  correlation  ( 4>ij ) 
and  total  turbulent  diffusion  ( T>ij )  of  Reynolds  stress: 

Gij  =  -p[(iH0)gj  +  (uj6)gi\  (2) 

Eij  —  2  v(UitkUjtk) 

<t>ij  —  (~~(ui,j  +  uj,i)) 

PO 

T)  T) 

'Dij  =  [— ( — Ui)5jk  —  (  uj)5ik  —  (uiUjUh)  +  i '  ('WiW-j  ),&],& • 

Po  Po 

The  acceleration  due  to  gravity  is  given  by  ft,  v  is  kinematic  viscosity  of  the  fluid,  and  po  is  the 
reference  density  of  the  fluid.  The  buoyant  production  and  dissipation  rate  of  turbulent  kinetic  energy 
(K  =  \{uiUi))  are,  respectively,  G  =  \Gu  and  e  =  \ea.  The  dissipation  rate  tensor  can  be  split  into 
its  isotropic  and  anisotropic  deviatoric  parts  as  follows: 

2 

Eij  —  —  +  dij.  (3) 

Citing  small  scale  isotropy,  the  anisotropy  of  dissipation  is  generally  neglected.  The  pressure-strain 
correlation  is  usually  modeled  as  [6] 

4 Hj  =  -Ci£(^^-  -  ^Sij)  +  C5f3((uie)gj  +  {u36)gt  -  |( uk6)gk5ij ), 


(4) 


where  C\  (usually  chosen  value  lies  between  3  and  5)  and  C5  (in  the  range  0.33  to  0.6)  are  numerical 
constants.  In  shear-generated  turbulence,  the  pressure-strain  correlation  model  will  also  be  a  function 
of  the  mean  flow  strain  rate  and  vorticity  tensors. 

The  evolution  equation  of  the  thermal  flux  is 

^  +  Uk(ui0)tk  =  Pie  +  Gie  +  (pie  +  P%e  —  (5) 

The  terms  in  the  equation  correspond  to:  the  time  rate  of  change,  advection,  thermal  production  f^), 
buoyancy  production  (Gie),  pressure-temperature-gradient  correlation  (pie),  total  turbulent  diffusion 
(Vie)  and  viscous  dissipation  (e^): 


Pie  =  -(uiUk)Tk;  Gie  = -Pgi(e2y,  (6) 

nr) 

(pie  =  (— 0,*);  £*0  =  (z'  +  k)  («;,*#,  fc) 

Po 

Vie  =  [-(—0)&ik  -  (uiOuk)  +  v(0uitk)  +  K(uid,k)],k- 

Po 

The  thermal  diffusivity  of  the  fluid  is  k.  The  pressure-temperature  correlation  (in  the  absence  of  mean 
flow  velocity  gradients)  is  usually  modeled  as  [6] 


4>ie  =  -Cle^(uid)+Czem(e2)). 


(7) 


The  value  of  the  constant  C\e  in  most  models  range  from  3  to  5  and  Cze  is  between  0.33  and  0.5. 
The  viscous  dissipation  (ei0)  of  thermal  flux  also  needs  to  be  modeled.  In  sufficiently  high  Reynolds 
number  flows,  the  small  scales  are  approximately  isotropic  leading  to  £io  being  zero.  At  lower  Reynolds 
numbers,  especially  near  the  walls,  viscous  dissipation  of  thermal  flux  may  be  non-zero  and  is  modeled 


as  [1] 


—  feeV^O- 


(‘ UiO ) 


(8) 


VkW)’ 

where  f£e  is  a  function  of  the  turbulent  Peclet  number  and  goes  to  zero  at  sufficiently  high  Reynolds 
number.  The  viscous  dissipation  rate  of  temperature  variance  is  given  bye#  =  n(0^,k)- 
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2.2  Algebraic  modeling  of  { UiUj )  and  (ut9) 


Often,  for  calculating  practical  flows,  it  is  computationally  too  expensive  to  solve  six  differential  equa¬ 
tions  for  the  Reynolds  stress  and  three  more  for  thermal  flux.  A  more  viable  approach  is  the  two- 
equation  level  of  turbulence  modeling.  At  this  level  of  turbulence  modeling,  the  mean  fields,  turbulent 
kinetic  energy,  scalar  variance  and  dissipation  rate  are  obtained  by  solving  their  respective  evolution 
equations.  Closure  is  achieved  by  modeling  Reynolds  stress  and  thermal  flux  with  algebraic  expressions 
in  terms  of  the  known  quantities.  We  now  seek  to  formulate  algebraic  models  for  the  anisotropy  of 
Reynolds  stress  and  the  correlation  coefficient  between  velocity  and  temperature  fluctuations.  We  will 
ultimately  verify  the  validity  of  the  modeling  assumptions  using  simulation  data. 

The  anisotropy  of  the  Reynolds  stress  is  defined  as 


b 


ij 


(uiUj)  1 

2  K  3 


(9) 


and, 

—^—1  =  2(Pq  —  £$)  +  Vq)  (13) 

at 

where  the  substantial  derivative  d/dt  =  d/dt  +  Uj{.),j.  In  the  above  equations,  the  production  of 
temperature  variance  is  given  by 

P$  =  —{iktyTi  (14) 
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and,  VK  =  \T>u  and  V0  =  (K(92)tkk  -  (uk02),k)  are  the  total  turbulent  transport  of  kinetic  energy  and 
temperature  variance  respectively. 

2.3  Transport  equations  for  bi3  and  Fi:. 


The  rate  of  change  of  the  anisotropy  tensor  and  thermal  flux  correlation  coefficient  can  be  evaluated 
using  chain  rule 


dbij  1  d(uiUj)  (uiUj)dK 

~dT  =  2K~~dt  2/if2  It  ’ 


(15) 


and, 

dFi  1  d(0Ui)  1.1  dK  1  d(92) 
dt  ~  JKW)  dt  2  i[K  dt  +  (i 9 2)  dt  J 


(16) 


Since  we  are  dealing  with  normalized  quantities,  it  is  best  to  consider  the  evolution  equations  in 
normalized  time.  The  turbulent  velocity  and  thermal  timescales  can  be  defined  respectively  as 


The  overall  timescale  (r) 
timescales: 


K  (92) 

Tv  —  5  Tq 

£  £0 


(17) 


of  the  buoyant  process  is  taken  to  be  the  geometric  mean  of  the  individual 


T  =  ^/TvT0. 


(18) 


Time  increment  in  normalized  time  is  given  by  dt*  =  dtjr. 

The  evolution  equation  of  the  anisotropy  tensor  can  be  easily  derived  by  substituting  equations  (1), 
(4)  and  (12)  into  equation  (15).  In  normalized  time,  the  evolution  equation  is 


dbij 

IF 


-\{Fig*  +  Fj9*  -  \Fkgt8ij)  -  \C^bl3 

4 — ^~Fi9j  +  FjQ*  —  — Fkgk5ij )  +  (Fkgk  + 

1  (nuj)  ,  r 

'■2KVij  ~  77772  Vk*  +  ~d 


2  K2 

§(!  -  y)  +  «K  + 


2K~‘lJ 
C5-  1 


(Fig*  +  Fjg*--Fk9t8ij) 


T\—Va 
[2K  13 


(lliUj)  t  ,  T 

-JkT-Pk)  +  w<hj. 


(19) 
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The  following  normalizations  have  been  used: 


rjr* 

1j 


K  T 


9*  =  P\ 


'm 

K 


rgi ■ 


(20) 


The  transport  equation  for  the  thermal  flux  correlation  coefficient  is  derived  by  substituting  equa¬ 
tions  (5),  (7),  (8),  (12)  and  (13)  into  equation  (16): 


f  =  +  ^ 


* 

V9i 


«(Ci*f  +  he)  -  '  1)  +  2t/l(§  -  1)1 


Tv 
Vie 


Tv"  e 

1  ppK  .  Ve  x, 

2  X  (02) 


Tg  £0 


lVkW) 

We  can  simplify  the  various  production  to  dissipation  ratios  as  follows: 


1  [rf'G 


1  pre-f  ~P9i(0ui), 


=  -29*f. 


Ty  r-Fg-i  _  /ly  _  _ P-'JP* 

tq  eg  V  Tfl  eg  J 


(21) 


(22) 


Weak  equilibrium  assumption.  In  slowly  evolving  flows,  the  mean  quantities  (T*  and  g*)  evolve 
more  slowly  in  space  and  time  than  the  turbulence  quantities  and  Ft) .  The  turbulence  quantities,  in 
normalized  time,  quickly  settle  down  to  equilibrium  values  on  the  imposition  of  mean  flow  parameters. 
The  weak-equilibrium  assumption  states  that  the  turbulence  is  approximately  in  equilibrium  with  the 
imposed-  mean  flow  parameters.  This  equilibrium  state  of  the  turbulence  can  be  obtained  by  setting 
the  substantial  derivative  of  the  turbulence  quantities  to  zero: 


dbij  _  dF{ 

dt*  dt* 


(23) 


Further,  in  the  weak-equilibrium  limit  we  assume  that  the  turbulent  transport  of  anisotropy  tensor 
and  the  thermal  flux  correlation  coefficient  are  negligible.  This  is  equivalent  to  the  following  models 
for  turbulent  transport: 

Vtj  =  ivK  (24) 

v»  =  + ^)\. 
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2.4  Algebraic  model  for  6^. 


The  transport  equation  (19)  for  the  anisotropy  of  Reynolds  stress,  on  invoking  the  weak-equilibrium 
assumptions  (equations  23  and  24)  yields 

-  g(C|-V-2F^(F^  +  F’S’‘  ~  (25) 

This  is  the  algebraic  model  for  the  anisotropy  of  Reynolds  stress  in  buoyancy  dominated  turbulent  flows. 
The  anisotropy  is  a  function  of  the  normalized  gravity  and  thermal  flux  vectors  and  the  constants  in 
the  pressure  strain-correlation  model. 

2.5  Algebraic  model  for  Fj. 


Invoking  the  approximations  in  equations  (23)  and  (24),  we  obtain  from  equation  (21) 


Fi  =  ~[2bijT;  +  ^T*  +  r]g*}, 


(26) 


where  77  =  1  —  C30  and 


Q  3  1[  J^(j  -  1)  +  -  *)]  +  Civ  +  1m- 

2  y  Ty  E  y  Tq  Eq  Ty 


(27) 


Using  equation  (22),  we  can  write 


Q  =  -Fi(T*  +  -g*)  +  Q0 


(28) 


where 

Qo  =  —  -)  —  i /-^  +  fee-  (29) 

V  Tv  2  V  Te 

The  expression  given  in  equation  (26)  is  not  yet  a  model  for  Ft  since  the  thermal  flux  appears  on  the 
right  hand  side  of  the  equation  also  (inQ).  The  modeling  will  be  complete  only  if  Q  can  be  expressed 
in  terms  of  known  quantities. 
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Determination  of  Q.  The  quantity  Q  must  be  such  that  the  model  for  F{  is  self-consistent.  Multi¬ 
plying  either  side  of  equation  (26)  by  (T*  +  ^g*)  we  obtain: 

Fi(T?  +  irf )  =  +  \t;  +  y>ijT’)(T-  +  ip*),  (30) 

which  can  be  restated  using  equation  (28)  as 

Q(Qo  -Q)  =  ~(T19i  +  \T*  +  2bijT*)(T*  +  l-gD,  (31) 

We  need  to  solve  this  equation  to  obtain  Q.  Such  a  value  of  Q  will  lead  to  a  self-consistent  model  for 
Fi. 

Define  the  following  invariants: 

h  =  «?(7?  + \glY,  h  =  irw  +  jSi)i  h  =  h,T-(T;  +  !«?).  (32) 

Substitution  of  (32)  into  equation  (31)  leads  to  a  quadratic  equation  for  Q: 

Q 2  -  QoQ  -  Jo  =  0,  (33) 

where 

h  =  vh  +  2^2  +  2/3-  (3^) 

Therefore, 

Q  =  i[Qo  ±  \/Qo  +  4/o]-  (35) 

Clearly  Q  has  to  be  real.  Therefore,  it  is  very  important  that  the  pressure-temperature  gradient 
correlation  model  coefficients  are  such  that  the  discriminant  is  always  positive. 

In  the  context  of  Rayleigh-Benard  convection  considered  here,  the  nondimensional  mean  thermal 
gradient,  T* ,  exists  only  along  the  vertical  direction  (T-j*  =  T2  =  0  and  T3  7^  0).  At  large  Rayleigh 
numbers,  and  correspondingly  large  Reynolds  numbers,  this  normalized  vertical  gradient  of  mean  tem¬ 
perature  is  large  and  negative  and  confined  to  thin  thermal  boundary  layers  near  the  top  and  bottom 
boundaries.  In  the  interior,  away  from  the  boundaries,  T3*  «  0.  The  nondimensional  gravity  vector  g\ 
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is  also  negative  since  it  points  opposite  to  the  vertical  z  direction.  The  negativity  of  both  T|  and  3* 
guarantees  that  in  Rayleigh-Bernard  convection  Jj,  h  and  I3  are  each  individually  positive  leading  to 
I0  being  positive  over  the  entire  convective  layer.  As  a  result,  the  discriminant,  Qq  +  4Jo,  is  always 
positive.  In  the  interior  Iq  can  be  approximated  as  f  (g*)2. 

Given  that  Iq  is  positive,  it  is  easy  to  show  that  one  of  the  roots  is  always  positive: 

Q(1)  =Qo  +  \jQl  +  4/o  >  0.  (36) 

The  other  root  is  always  negative: 

q(2)  =  Qq  —  yJ~Ql  +  4/o  <  0.  (37) 

On  inspection  of  the  model  for  thermal  flux  (equation  26)  it  is  clear  that  a  positive  value  of  Q  leads 
to  a  gradient  diffusion  of  thermal  flux,  whereas  a  negative  value  means  counter-gradient  diffusion. 
In  a  homogeneous  turbulent  convection  flow  near  the  ‘weak-equilibrium’  state,  counter-gradient  is 

unphysical.  Therefore,  we  deem  that  the  only  physically  permissible  model  for  Q  is  the  positive  root 
given  in  equation  (36).  Equation  (26)  in  conjunction  with  (36)  is  the  fully-explicit  and  self-consistent 

model  for  the  normalized  thermal  flux. 

3  Numerical  simulation  of  Rayleigh-Bernard  turbulence 

We  consider  the  classical  problem  of  Rayleigh-Benard  convection  in  a  layer  of  fluid  bounded  between 
two  horizontal  plates.  When  the  bottom  plate  is  maintained  sufficiently  hotter  than  the  top  plate, 
thermal  instability  drives  the  flow  and  at  large  enough  temperature  difference  the  flow  becomes  fully 
turbulent.  Since  in  this  configuration  the  mean  velocity  is  identically  zero,  turbulence  is  purely  driven 
by  buoyancy  and  shear  generation  is  absent.  A  DNS  database  ([7],  [8])  will  be  used  to  test  the  validity 
of  the  various  modeling  assumptions  and  the  effectiveness  of  algebraic  modeling  of  Reynolds  stress  and 
thermal  flux  terms  in  this  purely  buoyancy  driven  flow. 

Numerical  simulations  were  performed  in  a  box  of  a  square  platform  with  height  to  width  aspect 
ratio  of  2\/2.  A  schematic  of  the  computational  geometry  along  with  the  coordinates  is  shown  in  figure 
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1.  The  governing  Boussinesq  equation  along  with  the  incompressibility  condition  were  solved  using 
spectral  methods.  The  top  and  bottom  boundaries  were  considered  to  be  isothermal,  impermeable  and 
stress-free,  while  the  horizontal  directions  were  periodic.  One  of  the  nondimensional  parameters  in 
this  formulation,  the  Prandtl  number,  was  chosen  to  be  0.72  corresponding  to  that  of  air.  The  other 
nondimensional  parameter,  the  Rayleigh  number,  was  chosen  to  be  6.5  x  106.  This  Rayleigh  number  is 
nearly  four  orders  of  magnitude  greater  than  the  critical  Rayleigh  number  of  657  and  this  places  present 
simulation  in  the  hard  thermal  turbulence  regime,  according  to  the  classification  of  Castaing  et  al.  [9]. 
The  computations  were  well  resolved  with  a  uniform  grid  of  96  points  along  the  horizontal  directions 
and  97  points  along  the  vertical  direction.  Computational  data  was  collected  over  a  long  duration  of 
more  than  40  eddy  turn-over  times,  defined  as  H/K 1/2,  where  H  is  the  height  of  the  convecting  layer 
of  fluid.  All  of  the  data  presented  below  are  nondimensionalized  with  a  length  scale  of  H ,  velocity  scale 
of  -J RaPr /2Nu(k /H )  and  temperature  scale  of  ^/4Nu:i/RaPrAT,  where  «  is  the  thermal  diffusivity 
of  the  fluid,  AT  is  the  temperature  difference  between  the  top  and  the  bottom  boundaries.  The 
Nusselt  number,  Nu,  for  the  present  simulation  is  about  23.  The  above  proper  scaling  [10]  differs  from 
the  conventional  dilfusional  scaling  by  factors  \J RaPrj2Nu  and  \j4NuA  /  RaPr  in  the  velocity  and 
temperature  Beales.  The  diffusion  scaling  is  well  known  to  result  in  very  large  nondimensional  velocity 
at  large  Rayleigh  numbers  and  the  above  proper  scaling  will  reduce  all  nondimensional  quantities  to 
order  one. 

The  periodic  boundary  conditions  along  the  horizontal  directions  lead  to  translational  invariance 
and  statistical  homogeneity  along  the  x  and  y  directions.  The  presence  of  top  and  bottom  boundaries 
introduces  inhomogeneity  along  the  vertical  direction.  As  a  consequence,  single  point  turbulence  statis¬ 
tics  are  functions  of  only  the  vertical  direction.  Even  in  the  vertical  direction  approximate  homogeneity 
can  be  expected  in  the  interior  sufficiently  away  from  the  top  and  bottom  boundaries.  In  figure  2,  the 
mean  temperature,  <  6  >,  the  mean  square  temperature  fluctuation,  <  92  >  and  the  turbulent  kinetic 
energy,  K,  are  plotted  as  a  function  of  the  vertical  coordinate.  It  is  clear  that  the  rapid  variation  in 
the  mean  temperature  is  limited  to  the  two  thermal  boundary  layers  adjacent  to  the  top  and  bottom 
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boundaries  and  that  the  mean  temperature  is  nearly  uniform  in  the  interior  90%  of  the  convecting  layer. 
Owing  to  the  stress-free  boundary  conditions,  the  horizontal  components  of  velocity  are  non-zero  at 
the  top  and  bottom  boundaries,  resulting  in  relatively  large  kinetic  energy  at  z=0  and  1.  Both  <  62  > 
and  K  are  also  nearly  constant  over  the  interior  50%  of  the  layer. 

4  Evaluation  of  closure  models  and  modeling  assumptions 

In  this  Section,  we  use  simulation  data  to  evaluate  models  and  modeling  assumptions.  Specifically, 
we  investigate  the  following  three  important  modeling  issues:  (i)  pressure  correlation  models;  (ii)  the 
evolution  equation  budgets  of  Reynolds  stress  anisotropy  and  normalized  thermal  flux  to  evaluate  the 
validity  of  algebraic  stress  modeling  assumptions;  and,  (iii)  behavior  of  the  turbulent  transport  terms. 
Owing  to  homogeneity  along  the  horizontal  directions  (x-y  plane),  the  Reynolds  stress  anisotropy 
(M  and  thermal  flux  correlation  coefficient  (F*)  are  both  functions  of  only  the  vertical  (z)  direction. 
Furthermore,  invariance  between  x  and  y  directions,  along  with  the  constraint  bu  =  0,  leads  to  bn  = 
&22  =  ~\hz-  From  symmetry  arguments  the  other  off-diagonal  terms  of  the  Reynolds  stress  anisotropy 
tensor  are  zero.  Similarly,  mean  thermal  flux  exists  only  along  the  vertical  direction  and  the  horizontal 
components  of  the  thermal  flux  vector  are  zero.  Therefore,  in  the  following  sections  we  perform  the 
model  evaluation  on  the  fc33  and  F3  components  only.  Models  will  be  first  evaluated  with  standard 
commonly  used  model  constants.  Using  the  DNS,  optimal  values  of  model  coefficients  will  also  be 
evaluated.  Here  the  optimal  value  of  the  model  coefficient  is  that  which  minimizes  the  root  mean 
square  of  the  difference  between  the  data  and  the  model. 

4.1  Reynolds  stress  models 

Pressure-strain  correlation.  Away  from  the  walls,  the  pressure-strain  correlation  is  the  most  im¬ 
portant  turbulence  process  that  needs  closure  modeling.  The  fact  that  this  phenomenon  is  present  in 

the  simplest  of  homogeneous  flows  as  well  as  the  most  complex  has  made  this  one  of  the  most  studied 
turbulent  processes  in  mechanical  turbulence.  The  fluctuating  pressure  in  buoyancy-driven  turbulence 
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can  be  divided  into  two  parts,  each  governed  by  a  Poisson  equation: 


P%  =  -UijUjj  (38) 

P%  = 

where  ps  and  pf  are  the  slow  and  fast  (buoyancy)  components  of  pressure.  Along  the  same  lines  the 
pressure  strain  model  is  decomposed  into  the  slow  and  the  fast  pressure-strain  correlation  terms.  The 
slow  term  is  also  called  the  return  to  isotropy  term,  since  on  the  removal  of  all  turbulence  generating 
mechanism,  this  term  returns  the  turbulence  from  an  initial  anisotropic  state  to  isotropic  state. 

The  most  comprehensive  pressure-strain  correlation  model  in  buoyant  turbulence  is  given  by  Ris- 
torcelli  et  al  [11].  This  model  is  constructed  using  joint  realizability  constraints  on  Reynolds  stress  and 
turbulent  thermal  flux  and  has  the  right  behavior  in  rotating  flows.  The  Ristorcelli  model  is  however 
quite  complicated  in  structure  and  the  model  most  commonly  used  in  practical  applications  is  as  given 
in  equation  (4).  When  appropriately  normalized,  as  they  appear  in  the  bij  evolution  equation  (19),  the 
fast  and  the  slow  terms  of  the  model  are 

2^4  =  \c^Flg]^9*Fj~\Fkgl6tJ)  (39) 

T  s  „  TE 

2 K*ij  ~  Cl2Kij' 

where  C\  and  C5  are  numerical  constants.  The  more  general  pressure-strain  model  involves  additional 
constants  C2-C4  and  are  typically  used  with  production  terms  of  mechanical  origin. 

In  figure  3,  the  ^>33  component  as  calculated  from  DNS  data  is  presented  as  a  function  of  vertical 
height  z.  The  standard  model  with  a  commonly  chosen  value  of  0.5  for  C5  is  also  plotted.  Clearly, 
the  model  reproduces  the  DNS  rasults  very  poorly.  Even  the  qualitative  trend  is  not  captured  by  the 
model.  It  should  be  pointed  out  that  the  model  does  not  contain  any  near-wall  correction  terms  and, 
hence,  cannot  be  expected  to  do  well  close  to  the  wall.  But  in  the  center  region  of  the  flow,  the  model 
presents  a  nearly  flat  profile,  whereas  the  data  clearly  indicates  a  parabolic  profile  with  the  maximum 
at  the  center.  In  search  of  better  agreement,  we  hypothesize  an  extended  fast  pressure  model  of  a  type 
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sometimes  used  in  mechanical  turbulence: 

2^4  =  C*iicGbii  +  \C^J+9tFj  ~  l FkgtSij),  (40) 

where  G  is  the  buoyant  production  rate.  The  optimum  values  for  C*  and  C5  are  determined  to  be  —3.0 
and  0.54  respectively.  The  optimized  model  is  also  plotted  in  figure  3.  As  can  be  seen,  the  agreement 
away  from  the  walls  is  very  good.  The  negative  value  of  C*  is  significant,  for  it  implies  that  in  energetic 
turbulence  the  model  will  remain  realizable.  This  is  because,  in  energetic,  turbulence,  production  is 
positive  and  this  term  has  the  effect  of  bringing  anisotropy  back  towards  zero  which  is  indeed  the 
function  of  the  pressure  strain  correlation,  to  redistribute  the  kinetic  energy  equally  among  the  three 
components.  Had  the  coefficient  been  positive,  the  new  term  would  have  increased  the  anisotropy 
beyond  the  bounds  of  realizability  in  highly  energetic  turbulence. 

The  evaluation  of  the  slow-term  model  is  performed  in  figure  4.  The  DNS  data  indicates  that  the 
slow-term  is  very  large  at  the  walls  and  gets  progressively  smaller,  attaining  its  minimum  near  the 
center.  The  standard  model  captures  the  trend  qualitatively,  but  quantitative  agreement  leaves  a  lot  to 
be  desired.  When  an  optimized  value  is  used  for  Ci,  the  agreement  improves  slightly  near  the  walls  at 
the  expense  of  poorer  agreement  at  the  center.  An  expanded  model  which  included  a  non-linear  term 
in  anisotropy  was  tried  without  much  improvement  in  the  agreement.  This  leads  us  to  the  conclusion, 
that  the  rmnn>rical  coefficient  C\  should  perhaps  be  a  variable  depending  upon  the  local  state  of  the 
turbulence. 

Algebraic  model  verification  and  budget  of  fey.  The  budget  of  fey  evolution  equation  is  now 
investigated  to  evaluate  the  validity  of  the  algebraic  modeling  assumptions.  For  the  weak-equilibrium 
assumption  to  be  valid  we  should  have  negligible  rate  of  change  of  anisotropy  following  a  fluid  particle 
and  negligible  turbulent  transport  of  anisotropy  so  that  the  balance  is  between  production,  dissipation 
and  pressure-strain  correlation.  The  turbulence  under  consideration  here  is  statistically  stationary  and 
has  no  mean  velocity  field  and  hence  there  is  no  mean  flow  advection  of  anisotropy.  Therefore,  algebraic 
modeling  assumption  will  be  valid  if  the  turbulent  transport  of  btJ  is  negligible.  Referring  to  equation 
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(19),  the  assumptions  can  be  stated  as 


£>*>.  «  $ 
~  U,tj 


0;  4  =  -(^-4), 


(41) 


The  definitions  of  these  terms  can  be  easily  gathered  from  equation  (19): 
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In  figure  5,  4>\z,  -(P^-e 33),  2?^-  and  d \j  are  plotted.  In  order  for  the  algebraic  model  to  be  successful  it 
is  required  that  <633  s=>  -(P33  —£33)  and  the  remaining  terms  be  negligible.  This  is  clearly  not  the  case. 
The  turbulent  transport  of  anisotropy,  although  smaller  than  the  other  terms,  is  not  entirely  negligible 
even  in  the  center  of  the  channel.  This  term  is  large  enough  to  result  in  significant  differences  between 
the  pressure-strain  correlation  model  and  -(P33  -  £33)  .  Therefore,  we  conclude  that  the  algebraic 
Reynolds  stress  model  (equation  25)  may  not  be  appropriate  for  this  flow.  The  anisotropy  of  dissipation 
is  also  shown  in  the  figure  and  it  is  negligible. 


Turbulent  transport  of  6y.  The  modeling  assumption  that  the  turbulent  transport  of  fey  is  neg¬ 
ligible  implies  that  the  anisotropy  of  turbulent  transport  is  identical  to  that  of  the  Reynolds  stress 


itself: 


V° 
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From  figure  5,  it  is  clear  that  such  is  not  the  case  and,  hence,  we  will  now  take  a  closer  look  at 
turbulent  transport.  It  can  be  seen  from  equation  (2)  that  the  turbulent  transport  can  be  classified 
into  pressure  transport,  transport  through  triple  correlation  and  viscous  transport.  In  figure  6,  the  three 

components  of  V\z  calculated  from  the  DNS  data  are  shown  as  a  function  of  the  vertical  height.  While 
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the  viscous  transport  is  reasonably  small,  the  other  two  are  certainly  not  negligible.  The  pressure 

transport  is  particularly  large  near  the  walls  and  significant  even  near  the  center  and  needs  to  be 
carefully  accounted  for  in  modeling. 

4.2  Thermal  flux  models 


Pressure  temperature  correlations.  The  pressure  temperature-gradient  correlation  can  again  be 
decomposed  into  slow  and  fast  parts  which  are  modeled  commonly  as  follows: 

(44) 


Dol  and  Hanjalic  [2]  found  that  the  commonly  used  models  were  better  suited  for  simulating  the 
temperature  pressure-gradient  terms.  In  other  words,  they  propose 

<45> 


■»  CMsJf)) 

UXi 

Here  we  will  be  verify  their  observation  using  the  present  DNS  data. 

In  figure  7,  the  slow-term  model  is  tested  against  the  pressure  temperature-gradient  correlation  and 
the  temperature  pressure-gradient  correlation  evaluated  from  DNS  data.  Result  for  the  only  non-zero, 


z,  component  is  shown.  Clearly  the  model  reproduces  the  behavior  of  the  temperature  pressure-gradient 
data  very  well  throughout  the  domain  of  comparison  including  the  near  wall  regions.  The  optimum 
value  of  the  coefficient  is  found  to  be  Cye  =  2.3.  In  figure  8,  similar  comparison  is  performed  with 
the  fast  term  model.  In  the  interior  of  the  flow,  the  DNS  data  shows  that  the  correlations  are  nearly 
constant  at  their  respective  values.  The  fast-term  model  also  exhibits  nearly  flat  behavior.  Depending 
upon  the  value  of  the  model  coefficient  chosen,  either  of  the  correlation  is  reproduced  well.  The 
optimum  value  for  matching  the  pressure  temperature-gradient  correlation  is  C39  =  0.22  and  that  for 
temperature  pressure-gradient  correlation  is  C35  =  0.44.  The  near-wall  agreement  in  the  fast-term  case 
is  not  as  good  as  in  the  slow-term  case. 
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Dissipation  of  turbulent  thermal  flux.  The  model  for  dissipation  of  thermal  flux  is  given  in 

Section  2  and  is  repeated  here  for  convenience: 

£<e  =  fseV^eFh  (46) 


where  fEg  goes  to  zero  in  sufficiently  high  Reynolds  number  flows.  The  dissipation  of  thermal  flux 
should  be  identically  zero  in  high  Reynolds  number  turbulence  when  the  small  scales  are  statistically 
isotropic.  But  in  a  flow  such  as  the  present  one,  the  Reynolds  number  is  not  high  enough  for  £«  to 
vanish.  The  dissipation  e3e  calculated  from  DNS  data  is  presented  in  figure  9.  It  has  its  highest  value 
near  the  wall  and  is  fairly  low  in  the  center  of  the  flow.  Near  the  walls  the  turbulence  is  least  isotropic 
had  hence  the  value  of  this  dissipation  is  high.  Away  from  the  wall,  the  flow  is  more  isotropic  leading 
to  lower  values  of  £39.  The  value  of  the  model  coefficient  can  be  estimated  from  the  DNS  data  as: 


feB  = 


1  £3  6 


(47) 


F%  \Jz£e 

The  value  of  f£g  thus  calculated  is  plotted  in  figure  9.  The  coefficient  appears  to  be  a  fairly  strong 
function  of  z  close  to  the  walls,  but  is  nearly  a  constant  at  0.7  near  the  center  of  the  flow. 


Turbulent  transport  of  Fi.  The  turbulent  transport  of  Fi  can  be  inferred  from  equation  (21): 

'Di e 


F[g=r 


-  ei  ('Em  a-  —  m 
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The  standard  assumption  is  that  'Dig  can  be  modeled  as 


Vie-  2 (Ui9^  K  +  (02)} 


(49) 


so  that  D^y  can  be  neglected. 

The  validity  of  this  assumption  is  verified  in  figure  10,  where  V3g  and  +  ypj]  are  plot¬ 

ted  using  DNS  data.  First  we  plot  V:ie  as  defined  in  equation  (6)  which  is  consistent  with  with  the 
interpretation  of  <j>3g  as  the  pressure  temperature-gradient  correlation  (equation  44).  Also  plotted  in 
this  figure  is  V-iB  computed  without  the  pressure-temperature  correlation  (see  equation  6)  and  this 
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formulation  is  consistent  with  the  interpretation  of  foe  as  the  temperature  pressure-gradient  correla¬ 
tion  (equation  45).  The  model  appears  to  qualitatively  capture  the  features  well  in  either  case  but 
quantitative  agreement  is  only  reasonable.  The  overall  agreement  seems  somewhat  better  with  the 
inclusion  of  pressure-temperature  correlation  in  the  definition  of  X>39,  as  it  appears  in  equation  6. 

Budget  of  Fi  evolution  equation.  Since  the  turbulence  is  statistically  stationary,  the  evolution 
equation  of  Fi  can  be  rewritten  as  (see  equation  21) 

Pie  +  ei$  +  $£  +  =  0.  (50) 

The  production,  dissipation,  pressure-strain  redistribution  and  turbulent  transport  of  Fi  groupings 
can  be  surmised  from  equation  (21).  The  2  component  of  these  quantities  are  shown  in  figure  11. 
The  production,  dissipation  and  pressure  correlation  terms  are  clearly  larger  than  the  transport  term. 
Nonetheless,  the  transport  term  is  not  negligible,  especially  close  to  the  walls. 

Due  to  the  overall  reasonable  predictions  of  the  pressure-temperature  correlation  models  and  some¬ 
what  diminished  size  of  turbulent  transport  in  the  center  regions  of  the  flow,  the  algebraic  model  may 
be  more  appropriate  for  the  thermal  flux  than  for  the  Reynolds  stresses.  The  algebraic  model  for 
/•\  derived  in  Section  2  is  compared  against  DNS  data  in  figure  12.  The  model  displays  a  bi-modal 
behavior  not  seen  in  the  DNS  data.  The  maximum  disagreement  region  here  coincides  with  the  rather 
large  disagreement  region  in  figure  7  between  the  fast-pressure  model  and  corresponding  DNS  data. 
Despite  this  disagreement,  due  to  fundamental  validity  of  the  algebraic  assumption  in  the  case  of  Fi, 
we  believe  that  a  reasonable  algebraic  model  is  quite  possible  for  F>  The  search  for  a  better  algebraic 
model,  however,  should  start  with  the  development  of  better  pressure-gradient  temperature  correlation 
models. 

5  Conclusion 

The  various  turbulent  process  of  second  order  closure  modeling  interest  in  a  buoyancy  driven  turbulence 

are  closely  examined  using  direct  numerical  simulations  data  of  Rayleigh  Bernard  convection.  This  flow 
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was  selected  for  the  investigation  since  it  offers  a  unique  situation  where  buoyant  turbulence  can  be 
studied  without  the  complicating  influence  of  shear  turbulence. 

Our  study  demonstrates  that  the  commonly  used  pressure  strain  correlation  models  perform  some¬ 
what  poorly  even  in  the  interior  parts  of  the  flow  away  from  the  wall  effects.  While  the  prediction  of  the 
rapid  part  of  the  pressure  strain  correlation  can  be  improved  with  the  addition  of  an  extended  term, 
the  ability  to  predict  the  slow  pressure  strain  correlation  remains  a  challenge.  Turbulent  transport 
appears  to  play  an  important  part  in  the  evolution  of  Reynolds  stress  anisotropy,  and,  consequently, 
the  algebraic  stress  modeling  assumptions  are  not  well  satisfied  by  this  flow. 

On  the  thermal  flux  modeling  side,  our  study  reaffirms  the  observation  of  Dol  et  al  [2]  who  state 
that  the  current  models  capture  the  behavior  of  pressure-gradient  temperature  correlation  rather  than 
the  pressure  temperature-gradient  correlation.  This  would  obviate  the  current  practice  of  splitting 
the  pressure  correlation  term  into  homogeneous  and  inhomogeneous  parts.  The  data  shows  that  the 
thermal  flux  dissipation  is  not  negligible  and  that  the  model  coefficient  in  nearly  a  constant  away  from 
the  walls.  The  turbulent  transport  appears  to  play  a  somewhat  smaller,  but  still  significant  role  in 
the  evolution  of  thermal  flux.  Overall,  the  likelihood  of  a  reasonable  algebraic  model  appears  more 
promising  in  the  case  of  thermal  flux,  but  the  development  of  such  a  model  can  only  come  from  improved 
pressure-correlation  models. 

Fully- explicit  and  self-consistent  algebraic  models  for  the  Reynolds  stress  and  thermal  flux  have 
also  been  derived  from  their  respective  evolution  equations  using  the  weak-equilibrium  assumption. 
Although  not  quite  applicable  for  the  present  flow,  we  expect  these  models  to  be  adequate  for  realistic 
flows  at  higher  Reynolds  numbers  where  the  assumptions  made  are  likely  to  be  more  valid. 
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Figure  2:  Variation  in  the  mean  temperature,  <  6  >,  mean  square  temperature  fluctuation,  <  92  > 
and  turbulent  kinetic  energy,  K  in  the  vertical  direction. 
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Figure  3:  Pressure-strain  correlation  fast  term:  DNS  data  vs.  model.  Both  a  standard  model  with 
C\  =  0.0  and  C5  =  0.5  and  an  optimal  model  with  C\  =  -3.0  and  C5  =  0.54  are  compared  with  the 
DNS  data.  The  optimal  value  for  the  constants  are  evaluated  by  minimizing  the  difference  between 
the  model  and  the  DNS  data  over  the  interior  50%  of  the  layer. 
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Figure  4:  Pressure-strain  correlation  slow  term:  DNS  data  vs.  model.  Here  again  the  results  of  both 
the  standard  and  optimal  pressure-strain  model  are  compared  with  the  DNS  data. 
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Figure  7:  Pressure-temperature  correlation  for  the  slow  term.  DNS  data  vs.  model.  The  DNS  data  for 
both  the  pressure-temperature  gradient  and  the  temperature-pressure  gradient  terms  are  shown. 


Figure  8:  Pressure-temperature  correlation  for  the  fast  term.  DNS  data  vs.  model.  The  DNS  data  for 
both  the  pressure-temperature  gradient  and  the  temperature-pressure  gradient  terms  are  shown. 
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Figure  9:  A  plot  of  dissipation  of  turbulent  thermal  flux,  eie  against  2.  Also  shown  is  the  model 
constant,  fE$  evaluated  from  the  DNS  data. 
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Figure  10:  Turbulent  transport  of  FL. 
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Figure  11:  Budget  of  b\  evolution  equation. 
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Figure  12:  Algebraic  model  for  F3. 
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